arXiv:1507.01225vl [physics.flu-dyn] 5 Jul 2015 


Computational and physical consequences of interaction 
of closely located simultaneous hydraulic fractures 


Ewa Rejwer'^’*, Aleksandr Linkov'’’'^ 

°'Rzeszow University of Technology, 

Al. Powstancow Warszawy 12, 35-959 Rzeszow, Poland 
^Saint Petersburg State Polytechnical University, 
Polytechnicheskaya st., 29, St.Petersburg 195251, Russia 
^Institute for Problems in Mechanical Engineering, Russian Academy of Sciences, 
V.O., Bol’shoy Pr., 61, St.Petersburg 199178, Russia 


Abstract 

Strong interaction of closely located, nearly parallel liydranlic fractnres and 
its inflnence on their propagation are stndied. Both compntational and phys¬ 
ical aspects of the problem are considered. It is shown that from the compn¬ 
tational point of view, when a distance between cracks is small as compared 
with their sizes, the system becomes ill-conditioned and nnmerical resnlts 
deteriorate. The physical conseqnence of the interaction consists in decreas¬ 
ing of the crack opening and even greater decrease of condnctivity. Then the 
resistance to flnid flow grows what resnlts in the propagation of only those 
fractnres, the distance between which is large enongh. 

The research aims to snggests a means to overcome the compntational 
difficulty and to improve numerical simulation of hydraulic fractures in shales. 
Numerical experiments are carried out for a 2D problem by using the complex 
variable hypersingular boundary element method of higher order accuracy. 
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The condition number of the main matrix of a system, the opening of cracks, 
the stress intensity factors, effective conductivity and resistance of a cluster 
are calculated and analysed. The results imply that to avoid computational 
instability, it is possible to replace a cluster of cracks by a single properly 
located crack. The results also evidently conhrm that the propagation of 
distant individual fractures rather than of an array of closely located cracks 
is to be expected in practice, what agrees with observations. 

Keywords: hydraulic fractures, cracks interaction, hypersingular boundary 
element method, effective conductivity 


1. Introduction 


Hydraulic fracturing, when used to recover gas from low permeable shale 
structures, stimulates the flow of a fluid by generating the opening of a large 
number of micro cracks. Their growth and interaction with hydraulic frac¬ 
tures and pre-existing natural fractures affect the hydrofracture propagation 
and the efficiency of hydrofracturing treatment. The interaction has been 
studied by many authors, e^. N.R. Warpinski and L.W. Teufel fl. V.F. 


Koshelev and A. Ghassemi 


2| (these authors a 


so analysed in-situ stresses 


along the natural fractures); X. Weng et ah js], who attempted to model 
propagation of hydraulic fractures influenced by multiple natural fractures. 
However, because of complexity of the problem, there is still a need to un¬ 
derstand arising difficulties and to develop means to overcome them. 

It is well-known that when the ratio distance-to-crack size decreases, the 
interaction becomes stronger. When cracks are nearly parallel and the ratio 
is small, the opening of a fracture induces compressive stresses normal to the 
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drops and some of t 

n ^ 

shadowing, e.g. [4| 


fractnre snrface. This decreases the opening of neighboring, closely located 
fractnres. As a conseqnence, the propagation speed of some of the fractnres 

lem even stop. This phenomenon, called the stress 
, strongly affects the prodnctivity of mnltiple fractnres 
propagating from a wellbore. Thns, a rigorons analysis of the shadowing 
effect, tending to optimize the distance between snccessive treatments, is of 
importance for the petrolenm indnstry. It is to be taken into acconnt for 
proper modeling hydranlic fractnring. 

The paper aims to stndy interaction of closely located, parallel cracks 
nnder static conditions. We focus on the consequences of the interaction. In 
Section 2, the methods used for modeling 2D and 3D hydraulic fractures, 
are briefly reviewed. Section 3 contains the study of the interaction of mul¬ 
tiple parallel cracks, symmetrically or non-symmetrically located at various 
distances. At the beginning, some details of the method used are presented. 
Numerical experiments are carried out by employing the complex variable 
(CV) hypersingular boundary element method, based on the CV hypersin¬ 
gular boundary integral equations, tailored to account for the asymptotic 
behavior of helds in a close vicinityof singular points, such as tips of cracks, 
multi-wedge points, etc. (see e.g. |[9|-[n|). For the main matrix of a system, 
the condition number, which is a measure of nearness to computational in¬ 
stability, is calculated. It is shown, that when the distance between cracks 
decreases, or a number of cracks in a considered area grows, the condition 
number grows exponentially. This indicates ill-conditioning of the main ma¬ 
trix, which hnally causes deterioration of numerical results. Further subsec¬ 
tions contain the study of the stress shadowing effect, the openings of cracks 
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in a cluster and the values of the stress intensity factor (SIF) at the tips of 
cracks, for various distances between the cracks. Based on the calculated 
values of the openings, the effective conductivity of a cluster of cracks is 
analysed. The conclusions on the growth of the effective resistance to the 
flow of a fracturing fluid are made. Employing the obtained data, we con¬ 
clude on the possibility to replace a cluster of cracks by a properly located 
single crack. This serves to decrease the condition number of the matrix to 
the level excluding deterioration of numerical results. The numerical results 
provide also conclusions on the expected propagation of individual fractures. 
It is stated that the results of numerical experiments are consistent with the 
held observations. Section 4 contains the summary of the results obtained. 


2. Brief review of the methods commonly used for modeling hy¬ 
drofractures 

In this section, 2D and 3D fracture models developed for understanding 
hydraulic fracture process are briehy reviewed. Among them, those concern¬ 
ing with the analysis of the stress shadowing effect, are of special interest to 
our theme. 


2D fracture models and methods. The classica 


KGD model W S.A. Khristianovitch, Y.P. Zheltov 12|, J. Geertsma and 


2D models are the 


F. de Klerk 


13j | and the PKN model by T. Perkins, L. Kern |l^ and R. 


Nordgren ISj. The KGD model assumes the plane strain state in planes 


perpendicular to the fracture front, and it is used to study fractures on the 
initial stage of the propagation. In the PKN model, the plane-strain state is 
assumed in planes parallel to the fracture front. The both problems refer to 
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a homogeneous medium and they are reduced to solving in time a system of 
ID nonlinear differential (for the PKN model) or integro-differential (for the 
KGD model) equations in a single spatial coordinate. Presently the models 
account for a pumping regime and Carter’s leak-off. 

An approach for modeling a 2D fracture in an inhomogeneous medium 


with a net of crac 
by P.A. Cundall 


cs employs the distinct element method (DEM), proposed 


16| for solving rock mechanics problems. The DEM can 


simulate rocks subjected to static or dynamic loading. A region is divided 
by a discrete fracture network (DEN) that separates a hnite number of dis¬ 
crete blocks. Fluid leak-off is accounted for automatically as fluxes in cracks 
branching from the major fractures. 

3D fracture models. 3D models serve to study fracture nucleation and 
growth in a 3D space. Some of them account for multi-layered structure of 
rock formation. Specihcally, the pseudo 3D model (P3D, e.g. 1^) modihes 
the PKN model, adding height variation along the fracture length. Later 
on, the P3D model was adjusted to studying tight-gas formations in the 
form of so-called unconventional fracture model (UFM, 3|, 8|). The UFM 
serves to simulate fracture propagation, rock deformation, and fluid flow in 
the complex fracture network created during a hydrofracturing treatment. It 
combines a number of parallel and orthogonal PKN fractures. To account for 
the stress shadowing effect, it is supplied with corrective coefficients found 
from particular benchmark solutions. This also serves to roughly simulate 
the interaction of hydraulic fractures with pre-existing natural fractures. The 
model requires calibration of its parameters. 

An alternative way to overcome the limitations of P3D model consists in 
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Q, Q)- These 


mod- 


employing 3D variants of the DEM and DFN 
els combine geomechanical models and held observations. They are developed 
in mining and civil engineering and they are also used for coupled processes, 
including stimulation of fractured rock mass by huid injection, propagation 
of hydraulic fracture and its interaction with DFN. Still the DEM and DFN 
models are computationally expensive and they require calibration of param¬ 
eters characterizing elastic, huid and fracture properties. The calibration is 
performed by using available benchmark solutions of appropriate problems. 

The models discussed have been successfully employed for modeling frac¬ 
ture geometry and/or interactions between hydrofractures and pre-existing 
natural fractures. However, they do not allow to distinctly evaluate mutual 
inhuence of nearly parallel cracks, when the distance between them becomes 
notably less than the crack in-plane sizes. In view of signihcance of such 
conhgurations, we focus on their detailed investigation. 


3. Interaction of closely located fractures 

Strong interaction of nearly parallel closely located fractures concerns 
with both the computational and physical aspects of the problem. On one 
hand, there appear boundaries at small distances between them, so that 
after discretization the nodal points of meshes become very close; in limit, the 
boundaries coincide. Hence any numerical method, not specially designed for 
accounting for the asymptotic behavior of a solution in a thin strip between 
the boundary surfaces, unavoidably deteriorates when the thickness of the 
strip is too small. The situation is similar to that considered in the theories of 
thin plates and shells. Direct solution of exact equations becomes impossible 
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with decreasing thickness of a shell. Meanwhile the asymptotic analysis and 
proper using of the specihc geometry provide radical simplihcations and very 
accurate solutions of important engineering problems. Actually the strip 
between nearly parallel closely located cracks is in the state similar to that 
in a plate of the same thickness. Consequently it may be expected that some 
features of the plate deformation and bending will be reproduced in the strip 
between cracks. Meanwhile the degree of similarity and the edge effects 
are uncertain, and they are to be thoroughly examined by using a proper 
numerical method. The later should provide accurate and stable results up 
to the thickness, for which the asymptotic behavior of a solution becomes 
evident. 

Note that the thickness, for which asymptotics start to dominate, is not 
known in advance. This impels using of (i) a method of high accuracy and 
(ii) stability control of numerical results obtained. Below we meet these 
two requirements and conclude on possible simplihcations accounting for the 
asymptotic behavior of the solution. 

On the other hand, there arise physical effects specihc for the hydrofrac¬ 
ture problem. In particular, the shadowing ehect, mentioned in many papers 
on hydraulic fracture simulation, manifests itself in drastic diherence of the 
stresses in the strip between nearly parallel closely located cracks as com¬ 
pared with stresses near the surface of a single crack. It is just one of the 
physical ehects, generated by the crack interaction. It is closely related to 
other ehects, some of which are even of greater signihcance for the fracture 
propagation. The analysis below tends to reveal and to quantify them. 
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3.1. Brief description of the method used for calculations 


The calculations are carried out by using the complex variable hypersin¬ 
gular boundary element method (CV H-BEM) for the equation H 

, Af) 
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z ^ L 


( 1 ) 

where L is the total boundary of all cracks in a cluster, in which thermal 
and/or mechanical values may experience discontinuity; z E L is the complex 
coordinate of a field point; r is the complex coordinate of an integration 
point; the symbol A in front of a value, denotes its jump across the contour; 
(p is the potential (displacements); qn is the normal component of the flux 
(traction vector) at an element of L with the normal n] k = —kf, kf is the 
conductivity (shear modulus) of a matrix; {aZ is the angle between an 
element dz (dr) and the a:-axis; the index ’’plus” (’’minus”) refers to the side, 
with respect to which the normal n is outward (inward). 

The CV H-BIE ([T]) refers to anti-plane deformation and consequently 
Acf) is the discontinuity of the shear displacement through a crack. Still 
the equation keeps the general features of crack interaction and it serves 
us to simplify the analysis before studying more computationally involved 
plane-strain and 3D problems. (An extension to the latter problems on the 
basis of the results obtained is in progress). To keep track with the hydraulic 
fracturing we shall associate the displacement discontinuity with the opening. 

We solve the problem by employing the CV H-BEM of higher accuracy. 
To this end, higher-order approximation of the density functions and account¬ 
ing for singular behavior of the solution near crack tips are used. Specifically, 
the approximation includes CV polynomials of an arbitrary order and the 








power function with an arbitrary rational exponent (5: 


fji^) = 


(2) 


where c is the local complex coordinate of the end point of a boundary 
element; for a standard straight element j = 0,1,2,..., and for a stan¬ 
dard circular-arc element j = 0, ±1,±2,...; the exponent (3 accounts for 
the asymptotics of the density function near the singular points (e.g. tips 
of cracks). Thus, for an ordinary element /9 = 0. For a singular element, 
(3 = m/n G (—1,1) is a rational fracture, which in general is found by using 
a standard subroutine (see e.g. |2l|). 

For each boundary element, the global coordinates are transformed to its 
local coordinates so that an element becomes standard After using the 
approximation ([2]) in the local coordinates, the CV singular and hypersingu¬ 
lar, entering the CV H-BIE ([1]), for an element are 


S,{z) = 


fj{T)dT 


H,{z) = 


(r - zY 


(3) 


b b 

Without loss of generality, we assume b, c to be start and end point of 
the standard boundary element, respectively. When the standard element is 
singular, the c point at the tip of crack represents a singular point. For a 
standard straight element b = —1, c = 1; for a circular-arc element with the 
angle 26*o, b = exp{—i9o), c = exp{i9o). Below we use straight elements only. 

Evaluation of the influence coefficients (E]) over ordinary and singular 
elements is easily performed by employing the recurrent analytical formulae 
(see e.g (l0|) 


Sj+i = zSj{z) + aj, Hj+i = zHj{z) + Sj{z), j = 0,1,2,..., (4) 


9 







Vi = kSj{z)-aj-i], H^-i = j = 0, -1, -2,(5) 

z z 

c 

where Uj = f (t — cYr^dr are constants evaluated analytically. It should be 

h 

emphasized, that the formulae (jl]), (jS]) for the influence coefficients are exact 
and they notably decrease the time expense, as compared with numerical 
integration. 

The starting integrals Sq and Hq in (|1]) and ([5]) are: 
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[T — crar 


T — z 


T^dr 
T — t' 


Ho{z) = 


T 
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T^dr 

(7^’ 


b b—c b b—c 

where t = z — c. They are also evaluated analytically. Specihcally, for an 
ordinary element (/3 = 0 ) they are: 

1 1 


So{z) = In ^^ + 7 ri( 5 ^, Ho{z) = 


z c — z 


b-z ^ ' b 

where 5^ = 1 for 2 ; G [b, c], 5^ = 0 for 2 ; ^ [b, c], 2 ; = \z\ exp (iy). 

For a singular element (/? 7 ^ 0), it is sufficient to consider the case when 
/5 > 0. If /9 < 0 (|/9| < 1), we may write (3 = fdi — 1 with Pi = fd + 1 > 0. 
Then and we arrive at the integrals S'_i(z), H^pz) with positive 

/3i (0 < /9i < 1 ). Analytical formulae for the starting integrals 5 * 0 ( 2 :), Hq{z) 
with positive P = m/n {n > m), have been given in 23|. We present them 
in the simplihed forms: 

0 


= 


b—c ^ 






+ Tiit]p5t, ( 6 ) 
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where t = z — c, 6t = I ioi t E [6 —c, 0], <5* = 0 for t ^ [6 —c, 0], t = |t| exp {irj), 
ts = |t| " exp {i {rj + 2s7r) /n) is a root of the complex value t. The root — c 
is the principal one (s = 0). 

Numerical tests and experience, gained to the date, show that these forms 
of the CV BEM provide accurate and stable results (see e.g. Q 


n. M. 0. 


22j). The results of numerical experiments, presented in this paper, are ob¬ 
tained by employing the CV H-BIE ([2]) with the second order approximation 
of the density functions (displacement discontinuities and tractions) and a 
square root asymptotics near the tips of cracks {(3 = 0.5). 

3.2. Calculation of the condition number 

Note that when the number of cracks in a structure increases, the prob¬ 
ability of configurations with almost parallel, closely located cracks is high. 
When the distance between parallel cracks tends to zero, the main matrix 
contains rows with practically the same influence coefficients. This leads 
to ill-conditioning of the matrix and causes significant computational diffi¬ 
culties. Then, as mentioned, a numerical method not specially designed to 
account for this geometrical feature unavoidably fails and numerical results 
deteriorate. Therefore, it is essential to trace the condition number (see e.g. 


24|). which is a measure of nearness to computational instability. 

As an example, we consider the clusters of V = 2, 3, 4, 6 and 8 straight 
parallel cracks of the same length L = 2/ = 0.15 for symmetrical cracks, 
and L = 2/ = 0.1 for non-symmetrical cracks, located in a square region 
with the side length equal to 1.0. The square is subjected to the anti-plane 
shear. The distance between the external cracks in a cluster is d (see Fig.l). 
In the numerical experiments, the distance d decreases from 0.6Z to 0.006/. 
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Figure 1: The schemes for TV = 2, 4 and 6 of a) symmetrically and b) non-symmetrically 
located, parallel cracks of the length L and the distance between external cracks d. 

The cracks are located either symmetrically (Fig.la) or non-symmetrically 
(Fig. lb) about the horizontal middle cross-section. Each crack is discretized 
by four, three node straight boundary elements: two ordinary elements near 
the center of a crack and two singular at the tips of cracks. In calculations 
we assume: zero tractions at crack surfaces; zero and unit displacements at 
the left and right vertical sides of the square, respectively; and zero shear 
traction = 0 at its horizontal sides. The condition number is calculated 
for different numbers N of cracks in a cluster and for various distances d. 

The condition number for symmetrically located cracks exceeds 10^°-10^^ 
when d/L is less than 0.04 (Fig.2). For non-symmetrically located cracks 
the interaction is weaker and the condition number grows slower (Fig.3), as 
compared with symmetrically located cracks. Still the condition number for 
them exceeds 10®-10^° when the ratio d/L becomes less than 0.04. 
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Figure 2: Condition number for symmetrically located cracks. 
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Figure 3: Condition number for non-symetrically located cracks. 
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It is essential, that in the both cases, when the distance d decreases and/or 
the number N of cracks in a cluster grows, the condition number grows 
exponentially. Moreover, based on the results of numerical tests, we could 
see that the destabilizing effect of the interaction between cracks becomes 
signihcant when the condition number exceeds the values 10^-10®. Then the 
computational error is well above the prescribed tolerance. 

These results refer to a single cluster of cracks. When modeling strongly 
inhomogeneous structures (i.e. shale structures, etc.) with a high crack den¬ 
sity, the number of clusters of closely located, nearly parallel cracks increases. 
According to the performed calculations, this may drastically complicate the 
numerical analysis of a problem. 

3.3. Analysis of stress shadowing, erack opening, stress intensity faetors and 
conduetivity 

The analysis of quantities characterizing mutual influence of the cracks in 
a cluster is performed for the same conhgurations as those considered in the 
previous subsection (Fig.l). Calculations show that in the strips between 
closely located cracks, the stresses on areas parallel to cracks surfaces are 
almost constant and they are the same as those on the surfaces themselves. 
This is the so called stress shadowing effeet. It results in almost uniform 
strains in the strips between the cracks so that mutual displacements of two 
surfaces of a strip is practically zero when the thickness of a strip is small 
(recall that for a uniform strain, these mutual displacements approximately 
equal to the product of the strain by the strip thickness). Hence it could be 
expected that the openings are localized only at the first and the last cracks 
of a cluster, while the openings of internal cracks are practically zero. The 
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calculations of the openings confirm this suggestion. 
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Figure 4: Crack openings along = 2, 4 and 6 symmetrically located cracks in a cluster. 


Fig.4 presents the distributions of the openings along each of the cracks 
in a cluster of 2, 4 and 6 parallel cracks. The cracks are located at the 
distance d = 0.01 {d/L = 0.067). Clearly, in accordance with the suggestion, 
the openings of internal cracks are much less than those of the external ones. 
Specihcally, the openings of the hrst and the iV-th crack approximately equal 
to wci = wc^ ~ ws/2, while the openings of the intermediate cracks tend 
to zero. The major input into the sum is that of the external cracks. 

This effect becomes stronger with growing number of cracks in a cluster 
with hxed distance d between the external cracks (see Fig.4). We conclude 
that when the distance between cracks decreases, the sum of the openings 
'^Ci tends to the opening of a single crack ws (for the considered ex- 
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ample, ws = 0.148 at the center of the crack). 

Note that the distances, at which the discussed asymptotic distributions 
of the stresses, stains and displacements in the thin strips between cracks 
become evident, are well beyond the range of computational instability of 
the method employed. The results start to deteriorate merely at the ratio 
d/L = 0.0003, 0.006 and 0.033 for iV = 2, 4 and 6, respectively. Meanwhile, 
in the examples the considered ratio was notably greater {d/L = 0.067). 

For comparison, the values of openings calculated for the cluster of iV = 2 
and 4 non-symmetrically located cracks, are presented in Fig.5 and Fig.6, re¬ 
spectively. We see, that at points ’’shadowed” by other cracks the openings 
are less. In particular, for a cluster of = 4 cracks the stress shadowing, in¬ 
duced by the internal cracks, strongly decreases the openings of the external 
cracks. 



X 


- 2_1 - 2_2 


Figure 5: Crack openings along N = 2 non-symmetrically located cracks in a cluster. 
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Figure 6: Crack openings along N = A non-symmetrically located cracks in a cluster. 

Nevertheless, in the case of non-symmetric cracks, the shadowing is weaker 
and consequently the results are more stable than in the case of symmet¬ 
ric conhgurations. The deterioration occurs when the ratio d/L is notably 
smaller. Specihcally, for = 2, 4 and 6 the results start to deteriorate when 
d/L = 0.0001, 0.00075 and 0.01, respectively. 

The distribution of the openings on a crack uniquely dehnes the SIFs at 
its tips. Hence, as clear from the results for openings, it can be expected that 
for the symmetrical conhgurations of Fig.l, the SIFs will be close to zero for 
the internal cracks, while for the external cracks they will be two-fold less 
than those for an isolated single crack of the same length. The results of 
Table 1 comply with this suggestion. 

They present the SIFs normalized by the value K/jj = cr°^\/^, cor¬ 
responding to the SIF at the tips of a single crack with the half-length 
I = L/2 = 0.075, located at an inhnite plate under unit anti-plane shear 
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Table 1: The values of the normalized SIF Km /Kjjj at the tips of iV = 2 and = 4 
symmetrically located cracks, for decreasing distance d between cracks. 


d/L 

N = 2 

iV = 4 

2_1 

2_2 

4_1 

4_2 

4_3 

4_4 

0.67 

0.8942 

0.8942 

0.7351 

0.5949 

0.5949 

0.7351 

0.33 

0.8384 

0.8384 

0.6839 

0.5357 

0.5357 

0.6839 

0.06 

0.7497 

0.7497 

0.5438 

0.2646 

0.2646 

0.5438 

0.03 

0.6664 

0.6664 

0.5357 

0.1372 

0.1372 

0.5357 

0.006 

0.5396 

0.5396 

0.5319 

0.0077 

0.0077 

0.5319 

0.003 

0.5229 

0.5229 

0.5214 

0.0015 

0.0015 

0.5214 

single crack 

1.014 

1.014 


stresses = 1-0. (Note that we consider a finite square region, for which 
in the case of a single crack = 1.014). It can be seen that for 

the symmetrically located cracks, the normalized SIFs at the first and iV-th 
crack tend to Kjjj/2, while the SIFs at each of the internal cracks tend to 
zero. Clearly the sum of the SIFs is close to the SIF of a single crack. Hence, 
the stresses at distances from the tips, exceeding the distance between the 
cracks, are practically the same as for a single crack located in the middle 
of a cluster. Still this equivalence does not mean that the conditions for the 
propagation of individual cracks are the same as well; they are quite different 
what has important implications discussed below. 

In the case of non-symmetrically located cracks, shown in Fig. lb, the SIFs 
at the tips, which are beyond the shadowed zones, tend to K^jj. These are 
external tips of cracks 2_1 and 2_2 for the cluster of two cracks, and these are 
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the external tips of the cracks 4_2 and 4_3 for the cluster of four cracks. The 
SIFs at the remaining (shadowed) tips tend to zero with decreasing distance 
between the cracks. 

The results for the SIFs imply that at the stages of the hydraulic fractur¬ 
ing, for which the fracture toughness is signihcant, simultaneous symmetric 
propagation of a number of closely located parallel fractures is practically im¬ 
possible for two reasons. Firstly, for such conhgurations, the driving SIFs are 
notably (at least two-fold) less than for a single crack. This makes symmetric 
propagation quite unfavorable. Secondly, if because of statistical fluctuations, 
one of the fractures overruns its closely located neighbors, it ’’shadows” the 
tips of the neighbors. Then SIFs at tips of the neighbors decrease to a level 
below the critical SIF, dehned by the rock strength. Such an effect is espe¬ 
cially signihcant at the stage of the fracture initiation what complies with 
the results for this stage obtained in the paper [25|. 

This conclusion holds for the toughness dominated regime of the fracture 
propagation. Meanwhile, for the viscosity dominated regime we need to take 
into account for the changes in the crack conductivity, rather than those in 
SIFs. 

In numerical simulations of hydraulic fractures, a fracturing huid is com¬ 
monly modelled as a Newtonian huid. (When a huid is non-Newtonian, it is 
assumed to be Newtonian with the dynamic viscosity /i dehned by the char¬ 
acteristic strain rate; when a huid contains proppant, the dynamic viscosity 
is prescribed from specially designed experiments as a function of the density 
of the slurry). Then the hux q is connected with the pressure gradient dp/dx 
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by the Poiseuille equation 



where k is the conductivity of a crack. It is essential that the conductivity 
strongly depends on the crack opening; for a single crack, it is proportional 
to the third degree of its opening ws'- k = kg = tc|/(12/r). For a system of 
N parallel cracks, the effective conductivity is the sum of individual conduc¬ 
tivities: 


N 



( 8 ) 


i=l 


where wci is the opening of the f-th fracture in a cluster. The resistance to 
a flow of a viscous fluid is inverse to k: r = 1/k. 

The results of the previous subsection serve us to evaluate the changes in 
the effective conductivity caused by the crack interaction. It is reasonable to 
compare the conductivity of a cluster kc with that of a single crack ks for 
various numbers N of the cracks in a cluster and various normalized distances 
d/L between the cracks: 



Fig.7 and Fig.8 present the results of calculations for cracks located in 
a cluster symmetrically and non-symmetrically, respectively. It can be seen 
that in all the cases the conductivity of a cluster is considerably less than 
that of a single crack. The ratio kc/ks decreases as 1/(A^^). 

The decrease of the conductivity is the consequence of the fact that in¬ 
dividual openings of closely located parallel cracks are notably less than the 
opening of a single crack. In view of the cubic powers w^. in (|8]), this leads 
to fast decreasing of the effective conductivity. Then the total resistance of a 
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Figure 7: The dependence between the conductivity of multiple, symmetrically located 
cracks and a single crack, for various ratio d/L. 
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Figure 8: The dependence between the conductivity of multiple, non-symmetrically located 
cracks and a single crack, for various ratio d/L. 
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cluster to viscous flow drastically increases as compared with the resistance 
of a single crack. As a result, the gradient of the pressure, needed for the 
transportation of a given volume of a fluid, is much less for a single crack of 
the width ws than the gradient needed to push the volume through a system 
of hner channels with the same sum of openings. Therefore localization 
of a prescribed flux in a single fracture appears to be much more preferable 
than its distribution in a number of hner channels with the same summary 
opening. 

We conclude that, for the viscosity dominated regime of fracture propa¬ 
gation, the propagation of a number of closely located cracks is physically 
unfavorable. It implies that only those fractures, which are sufficiently far 
from other fractures of similar sizes may propagate. This conclusion agrees 
with the results^f numerical simulations and held observations presented in 
the papers [^, 


4. Conclusions 

The conclusions of the study are summarized as follows. 

(i) Employing the CV H-BEM, with the second order approximation of 
the density functions and square root approximation of the opening near 
crack tips, may serve for studying closely located, nearly parallel cracks up 
to distances, at which the asymptotics for thin strips between the cracks 
become evident. 

(ii) When the distance between cracks decreases or/and the number of 
cracks in a cluster grows, the crack interaction unavoidably leads to deteri¬ 
oration of numerical results. For the CV H-BEM used the results start to 
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deteriorate when the condition number of the main matrix exceeds 10^ —10®. 

(iii) The stress shadowing effect in a cluster of closely located nearly 
parallel cracks appears in actually uniform stresses and strains in thin strips 
between cracks. 

(iv) The openings of internal cracks in a cluster are nearly zero, while 
the openings of two external cracks are close to one-half of the opening of a 
single crack of the same size. Thus the sum of the openings in a cluster tends 
to the opening of a single crack when the number of internal cracks grows. 

(v) The distribution of SIFs in a cluster is similar to that of the openings 
because the SIFs at tips of a crack are completely defined by the opening of 
the crack. Specifically, the SIFs of internal cracks are nearly zero, while the 
SIFs of two external cracks are close to one-half of the SIFs at the tips of a 
single crack of the same size. Thus the sum of the SIFs in a cluster tends to 
the SIF of a single crack when the number of internal cracks grows. 

(vi) The fields of stresses, strains and displacements around a cluster 
of closely located nearly parallel cracks are close to those around a single 
crack located in the middle of the cluster. Consequently, the interaction 
between clusters may be modelled by placing merely one crack at the middle 
of each cluster. This approach may serve to avoid computational instability 
commented in the point (ii). 

(vii) The specific distribution of the SIFs at tips of interacting cracks 
leads to propagation of a single fracture rather than to simultaneous prop¬ 
agation of a number of closely located fractures in the toughness dominated 
regime. Similar effect occurs in the viscosity dominated regime due to the 
specific distribution of the openings^ which results in drastic decrease of the 
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conductivity of a cluster as compared with a single crack. 

(viii) Further investigation and an extension to clusters in 3D are needed 
to make quantitative recommendations for proper design of multiple hy¬ 
draulic fractures in low-permeable shales. 
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